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Abstract 

Atherosclerosis, the leading death in the United State, is a disease in which a plaque builds up inside the arteries. As the 
plaque continues to grow, the shear force of the blood flow through the decreasing cross section of the lumen increases. 
This force may eventually cause rupture of the plaque, resulting in the formation of thrombus, and possibly heart attack. It 
has long been recognized that the formation of a plaque relates to the cholesterol concentration in the blood. For example, 
individuals with LDL above 190 mg/dL and HDL below 40 mg/dL are at high risk, while individuals with LDL below 100 mg/ 
dL and HDL above 50 mg/dL are at no risk. In this paper, we developed a mathematical model of the formation of a plaque, 
which includes the following key variables: LDL and HDL, free radicals and oxidized LDL, MMP and TIMP, cytockines: MCP-1, 
IFN-y, IL-12 and PDGF, and cells: macrophages, foam cells, T cells and smooth muscle cells. The model is given by a system 
of partial differential equations with in evolving plaque. Simulations of the model show how the combination of the 
concentrations of LDL and HDL in the blood determine whether a plaque will grow or disappear. More precisely, we create a 
map, showing the risk of plaque development for any pair of values (LDL,HDL). 
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Introduction 

Atherosclerosis, hardening of the arteries, is the leading cause of 
death in the United States, and worldwide. The disease triggers 
heart attack or stroke, with total annual death of 900,000 in the 
United States [1] and 13 million worldwide [2]. 

Atherosclerosis is a disease in which a plaque builds up inside 
the arteries. A plaque contains low density lipoprotein (LDL), 
macrophages, smooth muscle cells (SMCs), platelets, and debris. 
The plaque constricts the lumen of the blood vessel thereby 
increasing the shear force of blood flow [3,4]. As the plaque 
continues to grow, the increased shear force may cause rupture of 
the plaque, possibly resulting in the formation of thrombus (blood 
clot) [3,5], ischemic stroke, and heart attack [3—5]. 

The process of plaque development begins with a lesion in the 
endothelial layer, allowing LDL, to move from the blood into the 
intima and becoming oxidized LDL (ox-LDL) by free radicals 
(FRs). FRs are oxidative agents continuously released by bio- 
chemical reactions within the body, including the intima [6-8]. 
Endothelial cells, sensing the presence of ox-LDL, secrete 
monocyte chemoattractant protein (MPC-1) [9,10], which triggers 
recruitment of monocytes into the intima [11]. After entering the 
intima, monocytes differentiate into macrophages, which have an 
affinity for the ox-LDL [12—14]. The ingestion of large amounts of 
ox-LDL transforms the fatty macrophages into foam cells [12,15]. 
Foam cells secrete chemokines which attract more macrophages 
[10,12,13]. SMCs from the media move into the intima by 
chemotactic forces due to MCP-1 [9,10], and platelet-derived 



growth factor (PDGF) [10,16], as well as by haptotaxis by the 
extracellular matrix (ECM). PDGF is secreted by macrophages, 
foam cells and SMCs [16,17]. ECM is remodeled by matrix 
metalloproteinase (MMP) produced by a variety of cell types 
including SMCs [18], and is inhibited by tissue inhibitor of 
metalloproteinase (TIMP) produced by macrophages and SMCs 
[19]. Interleukin IL-12, secreted by macrophages and foam cells 
[10,12,20], contribute to the growth of a plaque by activating T 
cells [9,20,21]. Indeed, the activated T cells secrete interferon 
IFN-y, which in turn activates macrophage in the intima 
[13,21,22]. At the same time that LDL enters the intima, high 
density lipoprotein (HDL) also enters into the intima, and becomes 
oxidized by free radicals [7,8] . However, oxidized HDL (ox-HDL) 
is not ingested by macrophages. HDL helps prevent atheroscle- 
rosis by removing cholesterol from foam cells, and by the limiting 
inflammatory processes that underline atherosclerosis [23]. Fur- 
thermore, HDL takes up free radicals that are otherwise available 
to LDL. Some of the key players in the atherosclerosis process are 
shown in Fig. 1. 

It has long been recognized that the cholesterol concentrations 
in the blood are indicators of the probability that a plaque will 
develop: higher LDL and lower HDL concentrations indicate a 
higher probability of plaque development. Public health guidelines 
in the U.S. specify what levels of LDL are low risk and what levels 
are high risk; they also specify what levels of HDL are poor and 
what levels are near ideal [24,25]. However, what is more relevant 
is to specify the risk associated with combined levels of LDL and 
HDL, and this is what the present paper addresses. A schematic of 
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Figure 1. Atherosclerosis schematics: the presence of ox-LDL in the intima causes monocytes to migrate from the lumen into the 
intima. Monocytes differentiate into macrophages which endocytose ox-LDL and become foam cells. SMCs are attracted from the media into intima 
by chemotaxis and haptotaxis. Cytokines released by macrophages, foam cells and SMCs activate T cells. T cells enhance activation of macrophages. 
HDL helps prevent atherosclerosis. 
doi:1 0.1 371 /journal.pone.0090497.g001 



the network of atherosclerosis is given in Fig. 2. In this paper, we 
developed a mathematical model of plaque formation by a system 
of partial differential equations based on Fig. 2. The aim of the 
model is to determine the risk of plaque formation for combined 
levels of LDL and HDL. In particular, we created a "risk-map" for 
plaque development in the LDL-HDL coordinate plane, where 
the first quadrant of the plane was divided into regions of high risk, 
low risk and no risk. Anti-cholesterol drugs are aimed at lowering 
high levels of LDL, but some drugs are known to also increase the 
level of HDL [24] . Hence such a risk-map may be important when 
evaluating the extend to which an anti-cholesterol drug can reduce 
the risk of atherosclerosis for particular individuals. 

Materials and Methods 

Mathematical model 

In this paper, we present a mathematical model based on the 
network shown in Fig. 2. The model includes the variables listed in 
Table 1. We assume that all cells are moving with a common 
velocity u; the velocity is the result of movement of macrophages, 
T cells and SMCs into the intima. We also assume that all species 
are diffusing with appropriate diffusion coefficients. The equation 
for each species of cells X has a form 

' l f+ViuX)-D x AX = F x , 
at 

where the expression on the left-hand side includes advection and 
diffusion, and Fx accounts for various growth factors, bio-chemical 
reactions, chemotaxis and haptotaxis. The equation for the 
chemical species are the same but without the advection term. 
Fig. 3 shows a 2D cross section of a blood vessel with plaque Q, 
and a planar cross section of a plaque in the direction along a 
blood vessel. 

Equations for lipoproteins [LDL (L), HDL (H), ox-LDL 
(io*)] and free radical (r). The distribution of LDL, HDL, ox- 



LDL and free radicals in the intima are described using reaction- 
diffusion equations [7], 

8 ^-D L AL= z kLrL, (1) 

reduction 



— -D H AH=^kHrH, (2) 

reduction 



—^-D Lox &L ox = kj,rL -), Lox mML ox , (3) 

production reduction 



- - D,Ar = r 0 - r(k L L + k H H), (4) 
ot v ' 

reduction 

where k L and k H are reaction rates of oxidization, and )il 0 xM is the 
reduction rate of ox-LDL due to ingestion by macrophages. Eqs. 
(1) and (2) model the evolution of LDL and HDL concentrations. 
It is assumed that LDL and HDL are lost by reaction of oxidation 
with free radicals. Equation (3) models the production of ox-LDL 
due to LDL oxidation by reaction with the radicals (first term on 
right-hand side) and a reduction of ox-LDL through ingestion by 
macrophages (second term on right-hand side). Equation (4) 
models the evolution of free radicals concentration with baseline 
growth r ( ). 

Equation for macrophages {M). The evolution of macro- 
phage density is modeled by 
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Figure 2. Schematic network of atherosclerosis. LDL and HDL are oxidized by free radicals, and become ox-LDL and ox-HDL respectively. Ox- 
LDL recruits macrophages to intima. By ingesting ox-LDL, macrophages are transformed to foam cells. SMCs are attracted into the intima by MCP-1 
(secreted by endothelial cells) and PDGF (secreted by macrophages and foam cells). Macrophages, foam cells and SMCs secrete IL-12, which activates 
T cells. IFN-y secreted by T cells enhance the activity of macrophages which contributes the plaque built-up. 
doi:1 0.1 371 /journal.pone.0090497.g002 



dM 



+ V-(uM)-D M AM-- 



-ViM Xc VP) + XuiyM 




(5) 



Here the first term on right-hand side accounts for recruitment of 
macrophages by MCP-1 [9], and the second term accounts for the 
activation of macrophages by IFN-y [13,21,22]. 

Equation for MCP-1 (P). The MCP-1 equation is given by 



— -D P AP= Aps-p — —j— 

St K Lox +L ox 



(6) 



production 



degradation 



where the first term on the right-hand side is the production of 
MCP-1 by endothelial cells, whose density is assumed to be 
constant, under the influence of ox-LDL [9] . 

Equation for T cells (J). The density of T cells, which are 
primarily CD4 + T cells [20], satisfies the equation 



8T 



+ S7(uT)-D T AT-- 




(7) 



In this equation, we assume that T cells are activated by IL- 1 2 in 
conjunction with MHC-II (major histocompatibility complex, class 
II). Actually, T cells are also activated by IL-1 and IL-6 produced 
by macrophages and SMCs [10,12,13]. However, because of lack 
of experimental data, we do not include the IL-1 and IL-6 
explicidy but instead consider their effect implicitly in estimating 



Table 1. The variables of the model: 


concentrations and densities are in units of g/cm 3 . 




L: concentration of LDL 


H: 


concentration of HDL 


L ox : concentration of ox-LDL 


r. 


concentration of free radicals 


P: concentration of MCP-1 


I v : 


concentration of IFN-y 


/ 12 : concentration of IL-12 


G: 


concentration of PDGF 


Q: concentration of MMPs 


Q r : 


concentration of TIMP 


M: density of macrophages 


T: 


density of T cells 


S: density of SMCs 


P- 


density of ECM 


F: density of foam cell 


(>: 


pressure (in g cm 2 /day) 


u: fluid velocity (in cm/day) 


doi:1 0.1 371 /journal.pone.0090497.t001 
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The first term of right-hand side is the production of I\% by 
macrophages enhanced by I y and resisted by HDL [23]. The 
production of I 12 by macrophages is resisted by 7 10 (which, for 
simplicity, is accounted by the factor Km 1 +m ) [10,12]. The second 

term represents the production of /12 by foam cells [20]. 

Equations for PDGF (G), MMP (Q) and TIMP (&.). We 

have the following sets of reaction diffusion equations for the 
chemokines (G, Q_and Qj): 



8G 
dt ' 



D G AG = X GM M + ). GF F + ). GS S_ ^doG 

production degradation 



(ii) 




M 



8Q 
dt 



D Q AQ-- 



production depletion degradation 



dt 



D Qr AQ r 



(12) 



a Qi . s S + a Qi . m M-- d Qr QQQ,. -dQ,.Qr^ (13) 

production depletion degradation 



Figure 3. Two 2D cross sections of a plaque. 1 M is the boundary of 
the intima in contact with the media, and T/ is the boundary of the 
intima in contact with the lumen. In (B) T L and r R are parts of the 
intima. 

doi:1 0.1 371 /journal.pone.0090497.g003 

the parameter ^r/ 12 - F° r simplicity, we include the anti- 
inflammatory effect of IL-10 produced by macrophages only 
implicitly, by the factor 1 /(Km + M). 

Equation for IFN-y (I ? ) . The dynamics of IFN-y concentra- 
tion is modeled by 



8I y 
It 



h y TT 



production degradation 



(8) 



where the first term on right-hand side represents production of Iy 
by T cells [21]. 

Equation for SMCs (S). The equation of the SMCs density is 
given by 



In equation (11), PDGF is produced by macrophages, foam cells, 
and SMCs [16,17]. In equation (12), MMP is secreted by SMCs 
[18] (first term on right-hand side), and is lost by binding with 
TIMP (second term). In equation (13), TIMP is produced by 
SMCs and macrophages [19]. 

Equation for foam cells (F). Macrophages that have 
ingested a large amount of ox-LDL become foam cells [7,12,15], 
so we have 



dF 
It" 



V(uF)-D F AF = X f 



M —dpF . (14) 



death 



production 



Equations for ECM (j>) and pressure (c). We assume that 
the intima has the constituency of a porous medium. Then, by 
Darcy's law, the velocity u of the cells is given by 



U= -Vf7, 



(15) 



r5 
77 



+ ViuS)-D s AS-- 



haptotaxis 



(9) 



The first two terms on right-hand side account for chemotaxis by 
MCP-1 [9,10], and PDGF [10,16], and the last term accounts for 
haptotaxis by ECM. 

Equation for IL-12 [In)- The concentration of IL-12 is 
modeled by 



M Iy ^ . F 

12 K M + M K In H + I, 11 K F + F^ 



(10) 



production 



degradation 



where a is the pressure. We also assume that the total density of all 
the cells plus the concentration of p is constant. This constant 
should be smaller than the average density of a plaque, 
1.22±0.03 g/ cm [26], because plaques contain some debris, 
which are not included in our model. We take the constant to be 
1 g/ cm , i.e., 



M+F+T+S + p = l. 



(16) 



We assume that all cells are approximately of the same volume 
and surface area, so that the diffusion coefficients of the all cells 
have the same coefficient, D. By adding Eqs. (5), (7), (9) and (14), 
we get 



dp 
dt 



-(l-p)Aa + V(7-Vp + DAp-- 



(17) 



where 
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<l,= -ViMx c VP) + lMi y M 
M 



" ~ h+K Iy 



-d M M 



+ A 



+ A 



12 Km + M 
V-(S Zc VG)-V-(5 to Vp) 



h2-d T T-V-(Sx c VP) 



(18) 



K 



-M — dpF. 



Eq. (17) gives a relation between p and a. We next derive an 
equation for p. The ECM is degraded by MMP [27], and is 
remodeled by macrophages and SMCs [18,27]. For simplicity, we 
take the remodeling rate to be a constant, X p , as in [28]. Then the 
equation of the density of ECM is given by 



8 -l +V-(up) = VC 1 " j)^qQp. 

degradation 

Since u = — V<7, this equation can be written in the form 



(19) 



dp 
dt 



-pAa-Va-Vp = tl/, 



where 



\l/ = A p p(l )-d pQ Qp. 

Po 



(20) 



Adding this equation and (17), we get our final equation for a: 

-Acr + Z)Ap = 0 + i//. (21) 
The equation for p can then be written as 



8 ^-DpAp + ((l> + il/)p-V(j-Vp = il/. (22) 



Boundary conditions 

For simplicity, we consider only 2-dimensional plaques as in 
Fig. 3. Then the boundary of the plaque consists of i) Fm, in 
contact with media; ii) a free boundary Fj, inside the lumen, and 
iii) two more vertical boundaries F L and T x of the intima in the 
case of Fig. 3(B). 

Boundary conditions on T/. We assume flux boundary 
conditions of the form 



except for M, and v.m(L 0X ,U) = <Xm j+TI' since ox-LDL attracts 
monocytes [11], while HDL limits the inflammation process[23]. 
Note that Lq and H a are the LDL and HDL concentrations in the 
blood, so we shall be interested to see how these concentrations 
determine whether a small plaque will grow or shrink. 

As in [29-31], we assume that the free boundary Fj is held 
together by cell-to-cell adhesion forces so that 

a = yK on Tj, 

where k is the mean curvature of the surface Fj. (If Tj is circular, 
then K is the reciprocal of the radius) Furthermore, the continuity 
condition un = V„, where n is the outward normal and V„ is the 
velocity of the free boundary Fj in the direction n, yields the 
relation 



Vn 



8a 



on Fr. 



(25) 



Boundary conditions on Y M . We assume non-flux boundary 
conditions for all variables except p and S on F M : For S, we have 



8S 
dn 



+ a. s (S-S 0 ) = 0, on F M 



where So is SMCs density in the media, and for simplicity we take 
tx s = a s (P,G) = a s ^^ since MCP-1 and PDGF attract SMCs 

from the media. As in the case of Fj, the boundary values of p are 
determined by Eq. (16). 

Boundary conditions on r^. and F R . We assume the 
periodic boundary conditions on F L and F R . 

Parameter estimation 

Table 2 lists the range of molecular weights of proteins and 
Table 3 lists their range of concentration. In the second columns in 
Tables 2 and 3, we indicate the (intermediate) values used in the 
simulations. The Tables 2 and 3 are used to estimate some of the 
model parameters. A summary of all the model parameters is 
given in Tables 4 and 5. 

Reaction rates. To estimate some of the parameters in the 
equations for proteins, we shall use the concept of "accessible 
surface area" [32,33] of a protein p, or briefly "area," A p , which is 
roughly the minimum surface area of the smooth shapes 
containing the protein. It was estimated in [34], that 
A LDL = l.02x 10~ n cm 2 , and A HDL = \3Ax 10~ 12 cm 2 , so that 
their ratio is 



^=4^=7.6. 



8X 



+ a x (X-X 0 ) = 0 



(23) 



for X= L, H, M, T, and non-flux boundary conditions for all other 
variables, 



8Y 

5iT : 



= 0, 



(24) 



where Y = L 0X , r, P, I y , S, I, 7 12 , G, Q Q, F. The boundary values 
for p are determined by Eq. (16). The coefficient Cf.x is a constant 



Accordingly, the corresponding reaction rates of the oxidation, k L 
and k H , are related by kp = 7.6kfj. Moreover, the reaction rate of 
oxidation of LDL by free radicals is kp = 10.37 gcm~ 3 day -1 
[7,34,35], so that k H = F36 g cm" 3 day" 1 . 

Diffusion coefficients. We assume that the diffusion coeffi- 
cients of all the cells are the same, and take them to be 
8.64 x 10~ 7 err? day 1 [28,36]. In order to estimate the diffusion 
coefficients of the various proteins, we assume that the diffusion 
coefficient of protein p, Dp, is proportional to its area Ap, i.e., 
Dp = KA p , where we take K to be the same for all small molecules. 
For glucose, which is a monomeric globular protein, A p can be 
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Table 2. Molecular weights. 





Protein 


Weight (kda) 


Explanation 


LDL 


549 


Over 95% of the LDL protein mass is apolipoprotein 


B-100 (apo B-100, 549 kDa (1000 g/mol)) [63]. 


HDL 


105 


The range of weight of HDL is 105-130 [64]. 


Free radical 


0.51 kda 


Free radicals include DPPH (0.39 kda), 


ABTS (0.51 kda) and superoxide anion (0.81 kda) [65]. 


IFN-y 


17 


IFN-y is described as a 17 kDa peptide [66]. 


PDGF 


35 


There are two PDGF polypeptides: 


PDGF-I with a molecular weight of about 35 kda, and 


PDGF-II with a molecular weight of about 32 kda [67]. 


MCP-1 


8.9 


[68] 


IL-12 


70 


[69] 


MMP 


52 


MMP-1 has two major species of molecular 


mass, 57 kDa and 52 kDa [48]. 


TIMP 


25 


The molecular weights of TIMP-1, TIMP-2 and TIMP-3 


are 28.5 kDa 21 kDa and 27 kDa respectively [48]. 


doi:1 0.1 371 /journal.pone.0090497.t002 



computed in terms of the molecular weight M p , by the formula 
A p = llAM^ 3 [33], and M ;; = 180 dalton [28], Z>,, = 1.04x 
10 _1 cm 2 day 1 [37]. Hence, for glucose, A" is determined by 

diffusion coefficient of glucose 



K-- 



Free radicals are monomeric globular proteins (average weight is 
500 da, Table 2). Hence 

D r = KA [Ke radicills =K x 1 1 . 1 X M\ '{I radicals 

= 2.93xl0 12 xll.lx(500) 2/3 xl0- 16 =2.05xl0- 1 cm 2 day-'. 

The diffusion coefficient of MMP is 4.32 x 10~ 2 cm 2 day" 1 [38]. 
We assume that the diffusion coefficient of TIMP is same as that of 
MMP. 

Production rates. We assume that, in Eq. (7), Ati u [I\t\~ 
2dx[T], where [X] denotes the average concentration of X We 
take [7i 2 ] «5 x 10- 10 gem" 3 , [7]*lx 10~ 3 gem" 3 from Table 3, 
compute D L = KA LDL = 2.93 x 10 12 x 1 .02 x anc j d T = 0.33 day" 1 [39,40] . Then 1 T , 12 is estimated by 1 x 10 6 

day . 



surface area of glucose 
diffusion coefficient of glucose 



n.iAf; 



2/3 



(26) 



1.04 x 10-' cm 2 day- 
11.1 x(180) 2/3 x 10- 16 



= 2.93 x 10 12 day' 



We can now 

10- 11 =29.89 cm 1 day" 1 and D H = KA HDL = 1>.9?> cm 1 day 



Table 3. Concentrations of proteins and cells. 





Proteins & cells 


Concentration (gem 3 ) 


Explanation 


LDL 


7x10~ 4 -1.9x10~ 3 


Range is 70-190 mg/dl [25,70]. 


HDL 


4x10~ 4 -6x10~ 4 


Range is 40-60 mg/dl [25,70]. 


IFN-y 


io- 9 


Range is of 0.1-10.0 ng/mL [10]. 


PDGF 


1.5x10~ 8 


Range in normal humans blood 


17.5±3.1 ng/mL [72]. 


MCP-1 


3x10- 10 


300 pg/ml [73] 


IL-12 


5x10" 10 


Range 200-800 pg/ml [20]. 


MMP 


3x10~ 8 


Range in plasma is 10—60 ng/ml [74]. 


TIMP 


3x10" 8 


Range in plasma is 10—60 ng/ml [74] 


SMC 


6x10~ 3 


Range 7,500,000-10,000,000 cells per ml [38]. 


Monocyte 


5x10" 5 


Range from 20,000 to 100,000 cells per ml [75]. 


T cell 


1 x10~ 3 


Range of CD4 + T cells in healthy normal adult 


of 500,000 to 1,500,000 cells per ml [49]. 



doi:1 0.1 371 /journal.pone.0090497.t003 
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Table 4. Parameters' description and value. 




Parameter 


Description 


Value 


k L 


reaction rate of LDL + Radical^>ox-LDL 


2.35x1(T 4 g~ 1 cm 3 day -1 [7,34,35] 


k H 


reaction rate of HDL + Radical— >ox-HDL 


5.29x1CT 6 g~' cm 3 day -1 [7,34] 


Dt 


diffusion coefficient of LDL 


29.89 cm 2 day -1 [33,34,37] & estimated 


D H 


diffusion coefficient of HDL 


3.93 cm 2 day -1 [33,34,37] & estimated 


Dl„ 


diffusion coefficient of oxidized LDL 


29.89 cm 2 day -1 [33,34,37] & estimated 


Or 


diffusion coefficient of radicals 


2.05x1(T 1 cm 2 day -1 [33,34,37] & estimated 


D M 


diffusion coefficient of macrophage 


8.64x1CT 7 cm 2 day -1 [28,36] 


D t 


diffusion coefficient of T-cell 


8.64x1CT 7 cm 2 day -1 [28,36] 


D, r 


diffusion coefficient of IFN-y 


1.08x10 2 cm 2 day -1 [76] 


As 


diffusion coefficient of SMCs 


8.64x1CT 7 cm 2 day -1 [28,36] 


D P 


diffusion coefficient of MCP-1 


17.28 cm 2 day -1 [44] 


Din 


diffusion coefficient of IL-12 


1.08x10 2 cm 2 day -1 [76] 


D G 


diffusion coefficient of PDGF 


8.64x1(T 2 cm 2 day -1 [42] 


Dq 


diffusion coefficient of MMP 


4.32x1(T 2 cm 2 day -1 [38] 


Dq, 


diffusion coefficient for TIMPs 


4.32x1(T 2 cm 2 day -1 [33,37,38] & estimated 


D f 


diffusion coefficient of foam cells 


8.64x1CT 7 cm 2 day -1 [28,36] 


^L ar M 


rate of ox-LDL ingestion by macrophages 


10 gem' 3 day -1 [7] 




activation rate of macrophages by IFN-y 


0.005 day -1 [39] & estimated 




production rate of MCP-1 


8.65x10~ 10 gem' 3 day -1 [44] & estimated 


*-TI n 


activation rate of T cells by IL-12 


1 x10 6 day -1 [39,40,49,74] & estimated 


h r T 


production rate of IFN-y by T cells 


0.066 day -1 [45,77] 


h xl M 


production rate of IL-12 by macrophages 


3x10~ 7 gem' 3 day -1 [45] 


hxiF 


production rate of IL-12 by foam cells 


1 x10~ 7 gem' 3 day -1 [45] & estimated 


^GM 


production rate of PDGF by macrophages 


0.1 day -1 [41,42] & estimated 


^GF 


production rate of PDGF by foam cells 


0.033 day -1 [41,42] & estimated 


"*-GS 


production rate of PDGF by SMCs 


0.5 day -1 [41,42] & estimated 


j-QS 


production rate of MMP by SMCs 


3x10~ 4 day -1 [28] 


A Q,S 


production rate of TIMP by SMCs 


3x10~ 5 day -1 [28] & estimated 




production rate of TIMP by macrophages 


6x10~ 5 day -1 [28] & estimated 




remodeling rate of ECM 


0.432 day -1 [28] 


^FM 


activation rate of foam cells 


0.12 day -1 [39] & estimated 


Am 


death rate of macrophage 


0.015 day -1 [39] 


dp 


degradation rate of MCP-1 


1.73 day -1 [44] 


dj 


death rate of T cell 


0.33 day -1 [39,40] 


*, 


degradation rate of IFN-y 


0.69 day -1 [37] 


ds 


death rate of SMC 


0.86 day -1 [7] 


dl„ 


degradation rate of IL-12 


1.188 day -1 [39,40] 


d G 


degradation rate of PDGF 


3.84 day -1 [42] 


d QQ, 


binding rate of MMP to TIMP 


4.98x10 s cm 3 g~' day -1 [44,48] & estimated 


d Q,Q 


binding rate of TIMP to MMP 


1.04x10 9 cm 3 g~' day -1 [44,48] & estimated 


da 


degradation rate of MMP 


4.32 day -1 [37] 


d Q , 


degradation rate of TIMP 


21.6 day -1 [46] & estimated 


d,, e 


degradation rate of ECM due to MMP 


2.59x10 7 cm 3 g~' day -1 [36] 


d F 


death rate of foam cell 


0.03 day -1 [39] & estimated 


doi:1 0.1 371 /journal.pone.0090497.t004 



PDGF is produced by SMCs, and likely also by endothelial cells i.e., A<jm = 0.1 day . Since SMCs production rate of PDGF is 

and macrophages [41]. In wound healing, macrophages produce higher than that by macrophages [41], we take Ags = 0.5 day - . 
PDGF at rate of 5.76 day [42]. Since the plaque formation is a The production rate of MMP by tumor cells was estimated in 

much slower process, we take this rate X GM to be much smaller, [28] to be 6 x 10~ 3 day -1 . We assume that SMCs produce MMP 
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at a much lower rate, namely, ).qs = 3 x 10~ 4 day -1 . Since SMCs 
produce MMP to enable them move into the intima by haptotaxis, 
we assume that they produce TIMP at a lower rate than MMP, 



and take lg, 



s= whs 



= 3x10 day . As macrophages pro- 



duce most of the TIMP [43], we take the production rate of TIMP 
by macrophages to be Ag r M = 2/l.g,s = 6 x 10~ 5 day . 

We assume that the production rate of MCP- 1 by endothelial 
cells, IpE, is twice that of dpPo, where Prj is the concentration of 
MCP in the blood, which is equal to 1 x 10~ 9 gcm~ 3 [44]. We 
assume that dp = \cIm =0.03 day -1 [39], and that, in Eq. (14), 
Afm[M] x4d F [F] and [F] « [M], so that A FM =4d F = 0.\2 day" 1 . 

By [45], Xj u m = 3 x 10~ 7 day -1 . We assume that foam cells 
have lower production rates of I 12 and PDGF than macrophages, 
and take Xj u F and Igp to be one third of the values of 1/ 12 m and 
^GM, respectively, so that X[ uF = 1 x 10~ 7 gcm~ 3 day and 
aqp =0.033 day . 

Degradation rates. The degradation rate of MMP is dg = 4.32 
day" 1 [37]. Since TIMP has a short half life compared to MMP [46], 
we take its degradation rate to be Jg r = 5dg = 2 1 . 6 day" 1 . 

In [47] , the binding rate of MMP and TIMP is reported to be 
3 x 10 s M~ l s~ l , where lM = the mass per mole, and the 
molecular weights of MMP and TIMP are 52 kda and 25 kda, 
respectively [48]. Accordingly, we derive the binding rate per 
Molar per second (by same formula as in [44]), 



3.0 x 10 s 



2 QQr 



A , x 52000 x 1.66 x lO" 24 

1000 cm 3 

= 4.98 x 10 8 cm 3 g~ l day' 1 , 



24 x 3600 



day 



and 



dQ r Q = 



3.0 x 10 s 



x 25000 x 1.66 x 10- 24 gx 



1000 cm 3 

= 1.04x 10 9 cm 3 g- 1 day- 1 , 



24x3600 



day 



where Na is called the Avogadro number, and is the number of 
molecular per dm. N A =6.02x 10 23 , and 1.66 x 10~ 24 g is the 
mass of a proton for atomic mass unit. 

Other parameters. The range of macrophages in the blood 
is 2x 10 s — 10 4 gem- 3 [75]; we take M 0 = 5 x 10~ 3 gem- 3 . 
The range of T cells in the blood is 0.5 xl0~ 3 — 
1.5 x 10~ 3 gem- 3 [49]; we take To = 1 x 10~ 3 gem- 3 . The range 
of SMCs is 5 x 1 0 3 — 8 x 1 0 3 gem — 3 [38]; we take 
S 0 = 6x 10~ 3 gem- 3 . We assume that K Lox is half of L 0 in Eq. 
(6), and similarly, K M = Mf- in Eq. (7), and K F =^in Eq. (10). We 
assume that the influx of LDL and HDL into the intima is larger 
than the influx of macrophages and SMCs, and take 
«i = «i/ = l-0o«~', and a s = 0L M = 0.2 cm- 1 . The influx of T 
cells is assumed to be smaller than that of macrophages, and we 
take o.t = 0.05 cm- 1 . 

Numerical methods 

Finite element implementation. In order to illustrate our 
numerical method, we consider the following diffusion equation 
with Robin boundary conditions: 



dX 
dt 



+ y V-(uX) - DAX = f(X) in CI, 



8X_ 



+ a(X-X I ) = 0 on T,, 



(27) 



5n 



+ p(X-X M ) = 0on r M , 



where yW(uX) is an advection term, and either y = 1 or y = 0 (no 
advection), and u= — Vc Multiplying the differential equation by 
an arbitrary function v, and performing integration by parts using 
the boundary conditions, we get 



dX 
, dt 



vdx+y 



XVa Vvdx+D [ VX-Vvdx 
n in 



(28) 



= [ f(X) vdx -y\ Xv ^ ds - [ vp(X - X M )d 
)n Jr> on J r M 

+ [ vaiX-X^ds. 



This is an equivalent formulation of the system (27), which is 
better suited for simulation. 

Similarly, Eq. (21) for a has the equivalent form: 



VaVvdx-D 



Vp-Vvdx-- 



((j> + \l/)vdx- 



V d £ds, (29) 



for an arbitrary function v. 

Galerkin discretization. The standard Galerkin discretiza- 
tion method uses finite dimensional subspaces Vf, to approximate 
the solution X. Let {cu,}^] be a basis of V/,, where J\ r is the 
number of nodes within the triangulation K. Let XfieV/, denote the 
numerical approximation of X at time t = ndt, where dt is the time 
step, X£ is written as a sum 



7 = 1 



(30) 



for coefficient c" to be determined. If ^ is approximated by 
X n+l_ X n 

-, then (28) is equivalent to 



dt 



(X n+l -X")coidx+dt[y 



X"V<TV(Oidx+D 



= dt( 



f(X")widx-y 



f da 

X"wi — ds 
Jr> <5n 



VX n+l -Vwidx] 



(31) 



C0iP(X n -X M )ds+ [ co i <x(X n -X I )ds), 



i;dx 



\ X n+1 Widx+dtD \ WX n+1 -W(OiC 

Jo. )(l 

= X n co i dx-ydt X n Va-Wm i dx+dt f(X")ajjdx (32) 
Jo. Jn Jsi 

+ dt \ (Dj(a(X" - Xj) - yX" ^ )ds - dt \ m${X n - X M )ds. 
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Table 5. Parameters' description and value. 



Parameter 


Description 


Value 




Xc 


chemotactic sensitivity parameter 


8.64 x 10" 1 -1.73 X 10 4 cm 5 g" 1 day" 


[36,37] (10)* 


Xh 


haptotaxis parameter 


8.64 x 10" 2 ~ 8.64 x 10 4 cm 5 g" 1 day" 


[36,37] (10 2 )* 


Lo 


source/influx of LDL in blood 


7x 10" 4 -1.9x 10" 3 gem' 3 [25] 




Ho 


source/influx of HDL in blood 


4xl0- 4 -6xl0- 4 £c?f7- 3 [25] 




'"0 


source/influx of free radical into intima 


0.26 gem" 3 day" 1 [34] 




M 0 


source/influx of macrophages from blood 


5x 10" 5 gem' 1 [75] 




T 0 


source/influx of T cells into intima 


1 x 10" 3 gcm~* [49] 




So 


source/influx of SMCs into intima 


fi-x 10~ 3 vrm~^ T711 




Po 


PCtiA Hcincitu 
L^J vi utrl 1 j 1 Ly 


1 x iu gem L^oj 




p 

i o 


MCP - 1 concentration 


i i f\— 10 nil 
5 X IU |y ij 




Go 


PDGF concentration 


1.5 x 10" 8 [72] 




U-L 


influx rate of LDL into intima 


1 .0 cm" 1 estimated 




<y-H 


influx rate of HDL into intima 


1.0 cm" 1 estimated 




a M 


influx rate of macrophage into intima 


0.2 cm" 1 estimated 




(J.T 


influx rate of T cells into intima 


0.05 cm" 1 estimated 




as 


influx rate of of SMCs into intima 


0.2 cm" 1 estimated 




Klox 


ox-LDL saturation for production of MCP-1 


0.5 gem" 3 [64] & estimated 




Km 


macrophages saturation for activation of T cells 


2.5 xlO" 5 gem" 3 [75] & estimated 




and production of IL-12 


K F 


foam cells saturation for production of IL-12 


2.5 xlO" 5 gem" 3 [75] & estimated 






IFN-}' saturation for activation of macrophages 


1 x10 -11 gem" 3 [39] 






IFN-y saturation for production of IL-12 


7x10 -11 gem" 3 [45] 





* Values chosen in the simulation. 
doi:1 0.1 371 /journal.pone.0090497.t005 



where d" is the vector (d"), B is a matrix (By), By = 
^VcOj-VcOj-dx, and <?/ = D J n (Vp-V(Oj)dx + J n (<j> + \l/)aijdx — 
L Wj ( 4- ds. 

Outline of the procedure. Suppose the domain Q(Z) has 
polygonal boundaries F/i) and r 0 (/). Then we can cover the 
closure Q(t) of Q(t) by a regular triangulation T of triangles, i.e., 
£i(f) = UteT T wnere each Tis a closed triangle. The triangular 
mesh, which is a basic thing that Finite Elements requires, is 
generated by distmesh [50], which is a mesh generation tool 
implemented in MATLAB, and our algorithm is outlined in 
Algorithm SI. For the detailed implementations, such as: 
construct basis functions over the triangulation, assemble the 
stiffness matrix, etc, see references [51,52]. 

Results 

Numerical simulation is initialized by a small formed plaque, 
(see Figs. 4, 5, 6). Five combined levels of LDL and HDL (Lo and 
Hq) are tested for 300 days: 

(a) (Lq,Hq) = (190,40): a small plaque doubles in size at 300 
days; 

(b) (Lq,Hq) = (150,45): a small plaque increases approximately 
50% at 300 days; 

(c) (Lq,Hq) = (130,50): a small plaque remains small at 300 
days; 

(d) (LoMo) = (110,52): a small plaque decreases approximately 
30% at 300 days; 



Recalling (30), we can rewrite the system (32) as a linear system 
of equations 



Ac" 



(33) 



where c" is the vector of (cj), and the coefficient matrix A = (Ay) 
and the right-hand side b = (bj) axe defined by 



Ay= VwrVwjdx + dt WjWjdx, 
Jn Jo. 



and 



bj- 



+ t 



(c"wj)widx — ydt (c"cOj)Va- Vwjdx 



f ( (tf(Oj))a)idx + dt co,(a( ^ (ejeo,) - Xj) 
2 i Jr / i 



8a 

-yX n — )ds-dt 
on 



0,fi(Y^(c-L»i)-X M )d S 



Similarly, setting a" t = YljL \ d"a>j, where a" t is a numerical 
approximation of a at time ndt, Eq. (29) can be written as follows 

Bd n = e, (34) 
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(e) (Lq,Hq) = (70,60): a small plaque almost disappear at 300 
days. 

Fig. 4 shows the growth of the plaque in case (a), Fig. 5 shows 
the shrinkage of the plaque in case (c), and Fig. 6 shows almost no 
plaque in case (e). In Fig. 7, the weight of the plaque, the 
summation of total cells, namely, J n (M + F + S + T)dQ, is plotted 
for these five scenarios of combined levels of LDL and HDL. 
Similarly to Fig. 7, we show in supporting information files 
how the populations of macrophages, SMCs, foam cells and T 
cells, as well as the concentration of ox- LDL, IFN-y and IL-12, 
vary for different levels of LDL and HDL shown in Figs. SI, S2, 
S3, S4, S5, S6, S7. 

Fig. 8 shows a risk-map of plaque development. To create the 
risk-map, we divided the LDL axis by 121 equidistant points, i.e, 
L, = 70 + i (? = 0,'--,120), and divided the HDL axis by 21 
equidistant points, i.e., Hj = 40 +j (J = 0, ■ ■ ■ ,20). For each pair 
(Lj,Hj), we computed the weight of the plaque, W(Li,Hj) after 
100 days on the domain shown in Fig. 3 (B), and formed the risk 
matrix 




Figure 4. Simulations for the atherosclerosis model of 300 days 
after an initial plaque is formed with W o = 40 mg/dL and 
Lq = 190 mg/dL. (A: Cross sections of a blood vessel, B:Cross sections 
along the blood vessel). 
doi:1 0.1 371 /journal.pone.0090497.g004 
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Figure 5. Simulations for the atherosclerosis model of 300 days 
after an initial plaque is formed with W o = 50 mg/dL and 
Z. o = 130 mg/dL. (A: Cross section of a blood vessel; B: Cross section 
along the blood vessel). 
doi:10.1371/journal.pone.0090497.g005 



_ W(Lj,Hj)— Wo 
iJ ~ W 0 

where Wo is the initial weight of the plaque. The vertical axis on 
the right of Fig. 8 shows the legend of the percentage of plaque 
growth or shrinkage. Accordingly, we divided the LDL-HDL 
plane into three regions: region I predicts high risk of plaque 
development, region III predicts no risk, and the intermediate 
region II predicts low risk. 

Sensitivity analysis 

In order to support the robustness of the simulation results, we 
ran sensitivity analysis on parameters which appear in the 
differential equations and in the boundary conditions. The 
parameters chosen are those whose baseline was somewhat 
crudely estimated while at the same time they seem to play an 
important role in the development of the plaque. Specifically, we 
chose all the 15 production rate parameters from the third box of 
Table 4, all the 5 influx rate parameters from the third box of 
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0.2 0.4 0.6 0.8 



Time: 100 



Time: 200 





50 



100 150 200 

Time (day) 



250 



300 



Figure 7. Plaque weights for different levels of LDL and HDL. 

The units of H 0 and L 0 are mg/dL 
doi:10.1371/joumal.pone.0090497.g007 

The most significant positively correlated parameters are Lo 
and its influx rate (I-l ■ This is not surprising since LDL initiates the 
plaque formation. The most significant negatively correlated 
parameters are Ho and its influx rate ajj. Indeed, since HDL 
reduces the availability of free radicals, it plays an important 
negative role in plaque formation. 



Time: 300 




Figure 6. Simulations for the atherosclerosis model of 300 days 
after an initial plaque is formed with H o = 60 mg/dL and 
i 0 = 70 mg/dL. (A: Cross section of a blood vessel; B: Cross section 
along the blood vessel). 
doi:1 0.1 371 /journal.pone.0090497.g006 

Table 5, and Lq, Ho. We list all these parameters with their range, 
baseline and unit in Table 6. 

Following the sensitivity analysis method described in [53], we 
performed Latin hypercube sampling and generated 100 samples 
to calculate the partial rank correlation coefficients (PRCC) and p- 
values with respect to the weight of the plaque after 300 days. The 
PRCCs are shown in Fig. 9, and all the p-values (not shown here) 
are less than 0.01. A positive PRCC (i.e., positive correlation) 
means that an increase in the parameter value will increase the 
weight of the plaque while a negative PRCC (i.e., negative 
correlation) means increase in the parameter will decrease the 
weight of the plaque. We note that Xqs is positively correlated, as 
it should be. Indeed, if aqs is increased then MMP (£7j is increased 
(Eq. (12)) so that ECM (p) is decreased (Eq. (19)) and hence the 
plaque weight \(M + F + S+ T)dQ is increased (Eq. (16)). As 
another example, note that ^l ox m is negatively correlated. Indeed, 
if ^l ox m is increased then L ox is decreased (Eq. (3)), and o;m in the 
boundary condition will decrease, leading to smaller M, and then 
to smaller T and F. Similar explanation can be given to the other 
parameters. 



Discussion 

Atherosclerosis is a disease in which a plaque builds up inside an 
artery. The process of plaque formation begins when, as a result of 
a lesion in the artery, cholesterols LDL and HDL enter the intima, 
and LDL becomes oxidized by free radicals. Upon sensing ox- 
LDL, endothelial cells secrete MCP-1 which attracts monocytes 
from the blood. As monocytes enter to the intima, they 
differentiate into macrophages that ingest the ox-LDL and become 
foam cells. Foam cells attract more macrophages, followed by T 
cells from the blood, and SMCs from the media. HDL reduces the 
available free radicals, as well as inflammation within the evolving 
plaque, thus HDL acts to block plaque growth. 

Public health guidelines in the U.S. specify that LDL level of 
100-129 mg/dL is near ideal, 130-159 mg/dL is borderline high, 
and 160-189 mg/dL is very high, whereas concentration of HDL 
above 60 mg/ dL is best, and below 40 mg/dL for men or below 
50 mg/ dL for women is poor [25] . An important question is how 
to evaluate the risk of atherosclerosis for a pair of LDL and HDL 
taken together. This question is addressed in the present paper. 
We built a mathematical model of plaque development by a 
system of partial differential equations. The model includes two 
parameters: ho, the level of LDL in the blood, and Ho, the the 
level of HDL in the blood. 

The model can simulate the evolution of a small plaque for any 
pair of values of (LqMo). In Figs. 4, 5, 6, we simulated the plaque 
evolution over a period of 300 days. For example, one extreme 
case of Lg = 190 mg/dL, Hq =40 mg/dL, the plaque doubled 
after 300 days; in another extreme case of Lq = 70 mg/ dL, 
Hq = 60 mg/dL, the plaque disappeared after 300 days. We 
created a risk-map by taking sampling points of LDL and HDL 
values, and computing the weight of the plaque for each pair 
(Lo,Hq) after 100 days. The map shown in Fig. 8, indicates the 
percentage of plaque growth or shrinkage for any such pair. We 
accordingly divided the (LDL,HDL) quadrant into three regions: 
high risk, low risk, and non risk. 

The need to consider the ratio of LDL/HDL in predicting 
coronary heart disease was suggested in a case study by [54]. The 
American Heart Association considers the ratio (Tc/HDL)>5 to 
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70 80 90 100 110 120 130 140 150 160 170 180 190 

LDL 



Figure 8. Risk map for plaque development: Region I high risk; Region II low risk; Region III no risk. The five points i /-,.//, ) whose 
plaque's weight was simulated in Fig. 7 over a period of 300 days are indicated by "x", 
doi:1 0.1 371 /journal.pone.0090497.g008 



Table 6. Parameters chosen for sensitivity analysis. 



Parameter 


Range 


Baseline 


Unit 




[5,20] 


10 


gcm~ 3 day -1 




[0.002, 0.01] 


0.005 


day" 1 


Ape 


[4.32 x 10" 10 , 1.73 x 10" 9 ] 


8.65 x 10" 10 


gem -3 day -1 




[5x 10 5 , 2x 10"] 


1 x 10 6 


day" 1 




[0.033, 0.132] 


0.066 


day" 1 




[1.5 x 10- 7 , 6x 10- 7 ] 


3xl0- 7 


gem -3 day -1 


A InF 


[5x 10- 8 , 2x 10- 7 ] 


1 x 10- 7 


gc-m~ 3 day -1 


A-GM 


[0.05, 0.2] 


0.1 


day -1 


A GF 


[0.016, 0.066] 


0.033 


day -1 


A-GS 


[0.25, 1] 


0.5 


day -1 


/-QS 


[1.5 x 1(T 4 , 6x 10-"] 


3xl0" 4 


day -1 


AQrS 


[1.5 x 1(T 5 , 6x 10~ 5 ] 


3xl0~ 5 


day -1 


A QrM 


[3xlO- 5 , 1.2x10-"] 


6x 10~ 5 


day -1 




[0.266, 0.864] 


0.432 


day -1 


AFM 


[0.06, 0.24] 


0.12 


day -1 


y-t 


[0.5, 2.0] 


1.0 


cm~ l 


Cl-H 


[0.5, 2.0] 


1.0 


cm 


5iK 


[0.1, 0.4] 


0.2 


cm~ l 


a T 


[0.025, 0.1] 


0.05 


cm 


«S 


[0.1, 0.4] 


0.2 


cm~ l 


io 


[5x 10-", 2x 10" 3 ] 


1 x 10" 3 


gcm~ 3 


H 0 


[2x 10-", 8x 10~ 4 ] 


5x 10-" 


gem' 3 



doi:1 0.1 371 /journal.pone.0090497.t006 
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Figure 9. The PRCC of parameters for sensitivity analysis. 

doi:1 0.1 371 /journal.pone.0090497.g009 

indicate high risk of heart disease, and the ratio (7c/ HDL) <3.5 
to be risk free [55], where Tc denotes the total cholesterol, which is 
calculated by the formula 

Tc = LDL + HDL + ~ 7> (7> = triglycerides) . 

Table 7 shows the National Cholesterol Education Program 
(NCEP) guidelines associated with plaque buildup [56]. Accord- 
ingly, for the five points (a)-(e) in Results we have: 

(a) j^5Z > 5 for any value of Tr; 

(b) > 5 if Tr is above normal; 

(c) 3.5 < < 5 if Tr is not very high; 

(d) Tfot < 3 . 5 if Tr is normal; 

(e) jysz <2-5 if Tr is normal;. 

According to the NCEP guidelines, (a) and (b) should be in the 
high risk region; (c) in the low risk region; and (d), (e) in the no risk 



region, as indeed they are according placed in the risk map in 
Fig. 8. 

Some anti-cholesterol drugs, such as statins, lower LDL and at 
the same time also increase the HDL [24] . It is important to know 
which drugs can best achieve the desired risk-free balance between 
LDL and HDL, that is, bring the individual's (Lq,Hq) into the no 
risk (or low risk) region. By focusing not on just reducing Lq or on 
just increasing 77o, but on moving the combined (Lq,Hq) to the no 
risk (or low risk) region in the shortest medically feasible path, we 
believe one could choose a more personalized medicine from those 
currently available, which will reduce the risk of atherosclerosis 
with the lowest amount of doze, thereby also possibly reducing 
potential negative side effects. 

To illustrate this approach, we note that for some drugs the ratio 
of decrease in LDL to increase in HDL is already known. For 
example, this ratio is 1/3 for the new experimental drug 
Evacetrapib. Some anti-cholesterol drugs only decrease LDL (e.g. 
Colestid) while others only increase HDL (e.g. Lofibra). We can 
the represent effect of such drugs by unit vectors in the (Lq,Hq) 
plane: for example, Colestid , Lofibra f , and Evacetrapid \ . 
In Fig. 10. we consider three individuals, A, B, and C in the high 
risk region. In order to move them to the low risk region with the 
minimum amount of medication (side effects are ignored), the 
individual should choose the drug for which the line segment from 
the individual initial position to the low risk region is the shortest 
(W e assume that the amount of drug is proportional to the length 
of the line segment). Thus the best drug for A is the one that 
primarily increases HDL. Similarly, C will do better with a drug 
that primarily decrease LDL, and B should use a drug with 
appropriate ratio of decreasing LDL to increasing HDL. 

Some work has been done on antioxidant therapy for reducing 
the risk of atherosclerosis, but so far it has had limited success in 
preventing cardiovascular diseases [57-59]. A review of studies in 
which antioxidant gene therapy has been successfully used is given 
in [60]. Our model could account for antioxidative medication 
once we gain a good understanding of how such medication affects 
the source of free radicals, r 0 , in Eq. (4). 



Table 7. National Cholesterol Education Program guidelines. 





LDL Cholesterol Level 


Category 


Less than 100 mg/dL 


Optimal 


100 to 129 mg/dL 


Near or above optimal 


130 to 159 mg/dL 


Borderline high 


160 to 189 mg/dL 


High 


190 mg/dL and above 


Very high 


HDL Cholesterol Level 


Category 


Less than 40 mg/dL (for men) 


Low HDL cholesterol. 


Less than 50 mg/dL (for women) 


A major risk factor for heart disease. 


60 mg/dL and above 


High HDL cholesterol. 


Triglyceride Level 


Category 


Less than 100 mg/dL 


Optimal 


Less than 150 mg/dL 


Normal 


150-199 mg/dL 


Borderline high 


200-499 mg/dL 


High 


500 mg/dL and above 


Very high 



doi:1 0.1 371 /journal.pone.0090497.t007 
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Figure 10. Drug treatment recommended for individuals A, B and C. 

doi:1 0.1 371 /journal.pone.0090497.g010 



Some of the parameters in the differential equations in our 
model had to be rather crudely estimated, since no data were 
available, while others may slightly vary depending on the 
individual. As more data become available, parameter values 
may be further refined. Our model uses only the values of LDL 
and HDL as biomarkers. It will be interesting in the future to 
incorporate also triglycerides into the risk-map. Future work 
should also explore how other risk factors, such as high blood 
pressure, smoking and diabetes affect the risk-map. 

We did not include in this paper the circulation ox-LDL in the 
blood, which is elevated only in patients with advanced 
atherosclerosis [61,62]. Our model could be extended to include 
this additional biomarker, but at present there is not enough data 
on how the level of ox-LDL in the blood correlates to a specific 
advanced state of the disease. 
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